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Abstract Using a Newtonian model of the Solar System with all 8 planets, we 
perform extensive tests on various symplectic integrators of high orders, searching 
for the best splitting scheme for long term studies in the Solar System. These com- 
parisons are made in Jacobi and Heliocentric coordinates and the implementation 
of the algorithms is fully detailed for practical use. We conclude that high order 
integrators should be privileged, with a preference for the new (10, 6, 4) method of 
( |Blanes etalj |2012[ ). 
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1 Introduction 



Due to their simplicity and stability properties, symplectic integrators have been 
widely used for long-term integrations of the Solar System, starting with the work 
of I Wisdom and HolmanM 19911). In many studies on the formation and evolution of 



the Solar System, where large numbers of particles are considered, the speed of the 
integrator is a major constraint and low order schemes have been often privileged 



as in the original scheme of (Wisdom and Holman 1991 ) or ( Kinoshita et all 1991 1 



(for a review see ( Morbidelli 20021). 
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On the opposite, in the present work we are focusing on high precision sym- 
plectic integrators that are designed for the computation of long term ephemerides 
of the Solar System, when one searches to reduce the numerical error of the algo- 
rithm to the level of the roundoff error of the machine. These integrators will also 
be useful for the detailed dynamical studies of the extra solar planetary system 
with strong planetary interactions. 

The first long term direct numerical integration of a realistic model of the 
Solar System, including all planets and the effects of general relativity and the 



Moon was made twenty years ago over 3 Myr (Quinn et al[ 19911 using a high 



order symmetric multistep method. This solution could be compared with success 
with the previous averaged solutions of ( |Laskar| |1989| 1990a ) and confirmed the 



existence of secular resonances in the Solar System ( Laskar et al 1992 1 . Soon 



after, using a symplectic integrator with mixed variables (Wisdom and Holman 



1991), Sussman and Wisdom| ( |1992[ ) could extend these computation to 100 Myr, 



confirming the chaotic behaviour of the Solar System discovered with the secular 
equations by Laskar (1989 1990a). 

As the Solar System is chaotic, the error in numerical integrations is multiplied 
by 10 every 10 Myr ( Laskar 1989 1 . Due to the limited accuracy of the models and 
initial conditions, it is thus hopeless to obtain a precise solution for the evolution 
of the Solar System over more than 100 Myr. The situation is even worse when 
the full Solar System is considered, as close encounters among the minor planets 
induce strong chaotic effects that will limit all possibilities of computing a precise 
solution for the planets to about 60 Myr ( Laskar et al| 2011a|b I. 

Despite this limitation, there is a strong need for precise ephemerides of the 
planets from the paleoclimate community. Indeed, the variations of the Earth or- 
bital elements induce some changes in the Earth climate that are reflected in the 
sedimentary records over million of years. This mechanism, known as Milankovitch 



theory (Milankovitch 19411 allows now to use the astronomical solution for the 



calibration of the geological time scales through the correlation of the variation 
of orbital and rotational elements of the Earth to geological records. This method 
has been successfully used for the Neogene period ( Lourens et al[ 2004) over 23 
Myr, and a large effort is pursued at present to extend this study over the full 
Cenozoic era, up to about 65 Myr. This quest led to search for high order sym- 
plectic schemes that are adapted to these long time computations, where high 
accuracy is requested |Laskar and Robutel (20011; Laskar et al (2004 2011a), but 
it should be noted that in the latest work, the integration of the Solar System 
model over 250 Myi[j] including five main asteroids took more than 18 months 
of CPU time. Some improvements of the algorithms were thus needed, and the 
present paper is the outcome of the studies that we have understaken in order to 
search for the best integrators for the next generations of numerical solutions. At 
the same time, we have compared various sets of coordinates (Heliocentrics and 
Jacobi), as the performances of these integrators depend on the choice of splitting 
of the Hamiltonian, and thus of the set of coordinates that correspond to these 
various splittings. As the integrators that are presented here are of high order, they 



1 Although it has been demonstrated that a precise solution o f the motion of the Earth 
cannot be computed over more than 60 Myr | |Laskar et al[[2011b^ , the solutions are system- 
atical ly computed over 250 Myr as some features of the solutions can be trusted over longer 
times |Laskar et al] ( |2004| |2011a| | . 



High precision Symplectic Integrators for the Solar System 



3 



can also be used for refined analysis of the newly discovered extra solar planetary 
systems, especially when the planetary interactions in the system are strong. 

For the planetary case, when using an appropriate set of coordinates, the equa- 
tions of motion are written as an integrable part Ha, that corresponds to the 
Keplerian motion of each planet, and a small perturbation Hg, given by the inter- 
action of the planets between each other. Hence, the system falls into the category 
of Hamiltonian system of the kind H = Ha + eHg. 

Several splitting integrating schemes that take advantage of this fact to de- 



1995 I was the first to 



lambers and Murison 



son | 
>12j| 



rive efficient integrators exist in the literature. McLachlan 
present such schemes and was followed independently by |C 
( |2000p and |Laskar and Robutel| ( |2001[ ). Most recently works by |Blanes et al| ( |2012 
derived higher order schemes that present very interesting behaviour. 

In this paper we describe these different splitting symplectic schemes and com- 
pare them for the case of the Solar System dynamics. We want to see which are 
the most efficient and accurate schemes. We will consider the gravitational N-body 
model and test the different integrating schemes against different planetary con- 
figurations, to be more specific: the 4 inner planets, the 4 outer planets and the 8 
planets in the Solar System (Section |4|. 

The Hamiltonian of the gravitational N-body problem H = T(p) + U(q) can be 
rewrite as H = Ha + eHg, using to different set of canonical coordinates: Jacobi 
and Heliocentric coordinates (Section [3]). The main difference between both set of 
equations is that in Jacobi coordinates the small perturbation Hg depends only in 
positions while in Heliocentric coordinates these one depends on both position and 
velocity. This is why in the literature Jacobi coordinates have been more widely 
used. In Section[5]we describe different symplectic schemes for Jacobi coordinates, 
and in Section [6]other symplectic schemes that are suitable for Heliocentric coor- 
dinates. In both Sections we describe and compare the different splitting schemes. 
Finally in Section[7]we compare the results for the two different set of coordinates. 



2 Splitting Symplectic Integrators (General Overview) 

Let H(q,p) be a Hamiltonian system where (q,p) are a set of canonical coordi- 
nates (i.e. q are the positions and p the momenta). It is well known that in many 
mechanical problems the Hamiltonian is of the form 

H(q,p)=T(p) + U(q), 

which is separable with respect to the local canonical coordinates. Using the Lie 
formalism we can write the equations of motion as: 

%={H,z} = L H z, (1) 

where by definition L x f := {x, f} is the differential operator L x , z = (q,p) and 
{■,•} denotes the Poisson brackelr] 

The formal solution of Eq. [1] at time t = tq + r that starts at time t = tq is 
given by 

z{r + t) = exp(rL ff )z(t ) = exp(r(L T + L v ))z{to). (2) 

2 rr r i BF 8G dF dG 

l r >"J — 2-ii=\ dpi B qi d qi dpi 
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In general the operators Lj< and L\j do not commute, exp(r(L^ + Ljj)) ^ 
exp(rLx) exp(ri[/), but we can find coefficients dj,&j such that for a given r, 



exp(r(L T + Lit)) = Y\ exp(airL T ) exp(6 i rL ;7 ) + 0(r r+1 ). 



(3) 



Using the Baker-Campbell-Hausdorff (BCH) identity we can find relations 
that the coefficients a^, bi must satisfy to have a high order scheme (Koseleff 1993 
1996). These are the so-called order conditions. For a given set of coefficients 
a,i,bi satisfying Eq. [3j the composition 



z{t) = S(t)z(t ) = Y\ exp(a;rL T ) exp(biTLu)z(To), 



(4) 



is a symplectic map of order r. 

The map 5(r) is symplectic because it is the product of elementary symplectic 
maps, exp(rLj') and exp(rL;y), and has order r because it approximates the exact 
solution up to order r r . We will refer to these kind of symplectic schemes as 
splitting symplectic integrators. 

Some of the main advantages of these kind of integrating schemes are: a) they 
are very easy to implement; b) they preserve the symplectic character of the Hamil- 
tonian system; and c) in general there is no systematic drift on the conservation 
of the energy during the numerical integration. 

These kind of symplectic schemes have been widely studied throughout the 



years by several authors (see Hairer et al (20061; McLachlan and Quispel (20021 
and references therein). As a matter of fact, splitting methods have been designed 
(often independently) and extensively used in fields as diverse as molecular dy- 
namics, simulations of storage rings in particle accelerators, quantum chemistry 
and, of course, celestial mechanics. 

There are several procedures to get the order conditions for the coefficients 
of the splitting scheme in Eq. [4j These are, generally speaking, large systems of 
polynomial equations in the coefficients that are obtained from Eq. [3] Two of the 
most popular are the recursive application of the BCH formula to the composition 
in Eq. |4j and a generalisation of the theory of rooted trees used in the analysis 
of RungeKutta methods due 



to 



Murua and Sanz-Serna ( 1999 ) (see also Hairer 
et al| ( p006[ )). The later procedure, while being more systematic than the former, 
is however not appropriate for the case splitting methods applied to Hamiltonians 
of the form A + eB. In Blanes et al (2012) a novel systematic way is proposed 



based on the so-called Lyndon multi-indices that is well adapted to that case. 

Splitting methods of order greater than two involve necessarily some negative 
coefficients and bj ( jGoldman and Kaper 1996 Sheng 1989 Suzuki 19911. 
Although this feature does not imply in principle any special impediment for the 
class of systems considered in this paper, it is clear that the presence of negative 
coefficients may affect the numerical error and the maximal step size of the scheme. 
For this reason, when dealing with high order methods, minimising the size of the 
negative coefficients and the sum of the absolute value of all the coefficients will 
be a critical factor in the choice of one particular set of coefficients. 

In this paper we do not intend to give the details on the derivation of the order 
conditions or how to find these coefficients. All these issues are analysed in detail 



High precision Symplectic Integrators for the Solar System 



5 



Blanes et al (2012). Our aim here is to compare the performance of different 



splitting symplectic schemes for the specific case of the integration of the Solar 
System. 

If we focus on the motion of the Solar System, or other planetary systems, we 
have a main massive body in the centre (the Sun) and the other bodies evolve 
around the centre mass following almost Keplerian orbits. We can take advantage 
of this to build more efficient schemes. Using an appropriate change of coordi- 
nates we can rewrite the Hamiltonian as, H = Hx + Hj (where \Hj\ <C \Hx\), a 
sum of the Keplerian motion of each planet around the central star and a small 
perturbation due to the interaction between the planets, where Hx and Hj are 
integrable. 



Wisdom and Holman (19911; Kinoshita et al (1991) where the first to split 



the Hamiltonian of the N-body problem in this way for numerical simulations of 
the Solar System, by means of what is called a mixed variable integrator, using 
elliptical coordinates to integrate the Keplerian motion. Splitting the Hamiltonian 
as Hx + Hj rather than the classical T(p) + U(q) already improves the performance 
of the leapfrog scheme. As \Hj\ is small with respect to l-ff/f |, the system falls into 
the class of Hamiltonian such that H = Ha + eHg for e small. In this particular 
case, the truncatio n order of the leapfrog scheme is no longer Ct 3 as for T(p)+U(q), 



but rather C'er 3 ( jMcLachlan 



1995 



Laskar and Robutel 



2001). 



In Sections [5] and [6] we will describe different families of symplectic splitting 
methods for Hamiltonian systems of the kind Ha + eHq and we will compare their 
performance for the particular case of the Solar System dynamics. 



3 The N-Body Problem 

Throughout this article we consider the non-relativistic gravitational N-body prob- 
lem as a test model for the different integrating schemes. We are aware that to 
have a realistic model for the Solar System dynamics one must include effects like 
general relativity or tidal dissipation. Nevertheless, and for the sake of simplicity, 
in this paper these effects are ignored as their presence should not compromise the 
performance of the schemes presented here. 

In a general framework, we consider the motion of n + 1 particles: the Sun 
and n planets, that are only affected by their mutual gravitational interaction. Let 
uo, ui, . . . , u n and uo, iii, • • • , u n be the position and velocities, in a barycentric 
reference frame, of the n + 1 bodies and let mo,mi,...,m„ be their respective 
masses. For simplicity, we consider mo to be the mass of the Sun and mi for 
i = 1, ... ,n the mass of the other planets. 

Taking the conjugated momenta U{ = m^ii;, the equations of motion are Hamil- 
tonian, with: 

H= 1 -y^- G y u mm \ r (5) 

2^ m; ^— ' Ui — Ui 

i=0 0<i<j<n J 

In this set of coordinates the Hamiltonian naturally splits into, H = T + U, where 
T depends only on the momenta (u ; ) and U depends only on the positions (u ; ). 

In general, when we deal with complex dynamical systems, it is important to 
take into account the relevant aspects of the system and use them to build efficient 
numerical tools to describe their dynamics. In the case of the Solar System we have 
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a massive body in the centre and the planets evolve following Keplerian orbits 
around it that vary through time due to their mutual interaction. 

Using an appropriate change of variables the Hamiltonian can be written as 
Hx + Hi, where \Hi\ is small with respect to |Hr-|, and both parts are integrable 
when we considered them on their own. There exist two canonical set of coordi- 
nates that allow us to split the Hamiltonian in this way: Jacobi and Heliocentric 
coordinates. 



3.1 Jacobi Coordinates 

The Jacobi set of coordinates have been widely used in Celestial Mechanics for 
developing analytical theories for the planetary motion. They where first used for 



the numerical integration of the Solar System by Wisdom and Holman ( 1991 1. 



Here the position of each planet, v; for i = 1, . . . ,n, is considered relative to 
the barycentre Gj_i of the previous i bodies, uo, . . . , Uj_i, and vo is taken as the 
centre of mass of the system: 



vo = (m u H + m n Un)/rin 

j=0 '■"j u j)/ r li~l 



v i = u i - (£-•-" m 



(6) 



where rji = J7J}=o m v have a canonical change of variables the momenta v; for 
i = 0, . . . ,n, must be: 



vo 
v; 



U H h G„ 



(7) 



In this set of coordinates the Hamiltonian in Eq. [EJtakes the form (Laskar 1990b I: 

1 Vi 



H 



Jb 



|vi|| 2 



2 Tjj_i m, 



G 



mjiii-i 
llvill 



+ G 



Vi-i 

|v; I 



l r il 



0<i<j<n 



(8) 



where Ai.j = ||uj — uj|| (the distance between the two bodies) can be expressed as 
a function of Vi and vj , and ri = Ui — uo . If we fix the centre of mass at the origin 
then vq = and vq = 0, and we reduce in 6 the number of equations of motion. 



3.2 Heliocentric Coordinates 

Here we consider the relative position of each planet with respect to the Sun: 



ro = u 

r ; = Ui - u 



and to have a canonical change of variables the momenta are: 



f o = u H h u n 1 

?i = Ui J 



(9) 



(10) 
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In this set of coordinates the Hamiltonian in Eq. [5] takes the form (Laskar 1990b I: 



mo + mi 
morrii 



(;■'"■■"■ ) ■ 

0<i<j<n 



G 



rn 



(11) 



where Aij = 1 1 r f — rj 1 1 for i, j > 0. If we consider the centre of mass of the system 
to be fixed at the origin we have that fo = 0, and we can easily recompute ro at 
all time. Hence, we have also reduced in 6 the number of equations of motion. 

One of the main differences between these two sets of coordinates is the size 
of the perturbation which in the case of Jacobi coordinates is smaller than for 
Heliocentric coordinates (see Table [l] in Section EL 

Moreover, in the case of Jacobi coordinates the perturbation part (Hj) depends 
only on positions so it is integrable when we consider it alone. But the expressions 
are more cumbersome than for Heliocentric coordinates (see Appendix [b]). While 
in the case of Heliocentric coordinates the perturbation part depends on both 
position and velocities, hence it is not integrable on its own. In Section [6] we will 
show how to adapt the splitting schemes to this particular case. 



4 Test to Perform 

Let S(t) = Yli=i exp(airA) exp(birB) be a splitting symplectic scheme. We say 
S(t) has s stages if it requires s evaluations of exp(rA) exp(r_B) per step-size. The 
smaller the step-size r used, the smaller is the error of the numerical solution pro- 
vided by the scheme, and the larger is the computational cost, as more evaluations 
of exp(rA) exp(r_B) are required to integrate over the same time period. 

Usually, the higher the order of the scheme the more number s of stages it 
requires, increasing the computational cost to advance a given step-size r. So a 
method with 4 stages will be more efficient than one with 2 stages if it can integrate 
a given accuracy with a step-size which is at least two times larger than the one 
required for the 2 stages scheme. In this sense, we define the inverse cost of S(r) 
as r/s, where s is the number of stages and r is the step-size used. Thus, if one 
scheme achieves the same precision than another scheme with smaller inverse cost, 
then we can say that the former is more efficient than the later. 

It is known that, for sufficiently small step-sizes r, the method S(r) integrates 
exactly (up to exponentially small errors that are below machine accuracy) a mod- 
ified Hamiltonian system that is close to the original one. Measuring the maximum 
variation of the energy along a given orbit will gives us a good idea of how close 
is that modified Hamiltonian to the original Hamiltonian. 

Motivated by that, in all our numerical test, we measure the relative precision 
of a given scheme applied with a given step-size r by computing the maximum 
variation (E^ = max{|jI(to) — H(t)\}) of the energy along a given numerical orbit 
obtained over 10 5 steps of the method (with the same initial conditions at the 
initial time to) and plot the Ei versus the inverse cost r/s. To fix criteria we will 
always consider step-sizes of the form: = 1/2 1 for i = 0, . . . , N . 

We are interested in very precise integrations of the Solar System, hence the 
main goal is to determine for each scheme the maximum step-size (r^) required to 
have an error in the energy variation up to machine accuracy. 
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Through the paper we consider three test models that we believe illustrate 
different particularities of the Solar System and can be extrapolated to other 
planetary systems. These are: a) the motion of the four inner planets (Mercury 
to Mars); b) the motion of the four outer planets (Jupiter to Neptune) and c) 
the motion of the eight planets on the Solar System (Mercury to Neptune). The 
initial conditions and mass parameters have been taken from the JPL Solar System 
ephemerides DEA405 (http://ssd.jpl.nasa.gov/). 

Table[l]shows estimates on the size of the perturbation for these three examples 
for both set of coordinates Jacobi and Heliocentric. To estimate the size of the 
perturbation we have integrated each system over 100 years and computed the 
maximum values for \Hj\ and \Hk\ along this integration. Here HKep represents 
the size of the Keplerian part and Hlmax the size of the perturbation part and the 
estimated size of the perturbation is given by s = Hlmax/HKep. 



Table 1 Size of the perturbation in Jacobi and Heliocentric coordinates for the three test 
examples considered in this work: 4 inner planets (Mercury to Mars), first line; 4 outer planets 
(Jupiter to Neptune), middle line; 8 planets on the Solar system (Mercury to Neptune), third 
line. 



Jacobi Coord. 
HKep Hlmax e 



1 . 3945E-04 
4.2924E-03 
4.4319E-03 



6.3342E-10 
8.7162E-07 
8.7158E-07 



4.5420E-06 
2.0306E-04 
1.9666E-04 



Heliocentric Coord. 
HKep Hlmax £ 



1.3945E-04 
4.2920E-03 
4.4314E-03 



9.1652E-10 
2.7184E-06 
2.8042E-06 



6.5720E-06 
6.3336E-04 
6.3281E-04 



We note that all the simulations in this article have been done using an ex- 
tended real arithmetics and that we use the compensated summation during the 
intermediate evaluation of exp(a^TA) and exp(fojT-B) (see Appendix |A|). 



5 Splitting Symplectic Integrators for Jacobi Coordinates 



In Section [3] we have seen that with an appropriate change of variables we can 
rewrite the Hamiltonian of the N-body planetary system as H^ + Hj where \Hj\ <C 
| -Hie | • Hence, the system falls into the class of Hamiltonian that can be expressed 
as 

H = H A + eH B , (12) 

with e <C 1. We can take advantage of this to build efficient high-order splitting 

In this sec- 



symplectic integrators (McLachlan 



1995 



Laskar and Robutel 



2001) 



tion we summarise the main ideas behind these methods and review some of the 
most relevant schemes. 



Using the Lie formalism the formal solution of Eq. ( 12 ) is 



!(t) = exp[r(A + eB)]z(t ), 



(13) 



where to simplify notation we use A = {Ha,-} = Ljj a , B = {Hg,-} = Ljj b - We 
recall that Ha and Hb are integrable, hence we can compute explicitly exp(r^l) 
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and exp(re_B). To have a splitting symplectic integrator of order r, we need to find 
the coefficients a,, 6, such that 

s 

S r {r) = Yl exp(a,rA) exp(e6. J rB), (14) 

8=1 

satisfies |5 r (r) - exp[r(A + eB)}\ = 0{r r+1 ). The Baker-Campbell-Hausdorff 
(BCH) theorem ensures us that S r (r) = exp(r"H), where H is also a Hamilto- 
nian system and belongs to the free Lie algebra generated by A and B, C(A,B). 
Moreover, we can express H as a double asymptotic series in r and e: 

tH= rp lfi A + erp hl B + et 2 P2 ,i[A, B] + ET 3 p 3A [[A, B], A] 

+ e 2 T 3 P 3,2[[A,B],B] + st a P4a[[[A,B],A),A] (15) 
+ e 2 T i P ^ 2 [[[A,BlB],A] + e 3 r A p 4 , 3 [[[A,B],B},B} + 

where ptj are polynomials in and b^. 

To have a symplectic scheme of order r we need: 

Pl,0 = at + a 2 + ■ ■ ■ + a s = 1, 
Pi,i = h + 62 + • • • + b s = 1, 
Pt, 3 = 0, Vi,j < r. 

The scheme 5 r (r) is symmetric if it verifies 5^ 1 (r) = S r {— r), in which case all 
the even terms in r in Eq. |15| vanish, having less conditions to satisfy for a scheme 
of a given order, r, enabling us to find high-order schemes at lower computational 
cost. There are two different types of symmetric compositions (Eq. 14): one in 
which the first and last exponentials correspond to the A part (and thus called 
ABA composition), 

ABA : e airA e £blTB e a2rA ■ ■ ■ e a2rA e eblTB e airA (16) 

and the other in which the role of exp(rA) and exp(er_B) is interchanged (BAB 
composition) : 

BAB : e sblTB e airA e £ll2TB ■ ■ ■ e £t>2TB e airA e eblTB . (17) 

All the integration schemes that we present in this paper correspond to the ABA 
class. For the experiments carried out, we have not found substantial differences 
in the efficiency with respect to methods in the BAB class. 

Notice that for symmetric methods, the last exponential at one step can be 
concatenated with the first one at the next integration step when the method is 
iterated, so the number of exponentials exp(rA) and exp(erB) per step is s, the 
number of stages. 

It is clear that in many cases \e\ <C r (or at least e « r). So we can have 
high-order schemes by only killing the error terms with small powers of e, and 
save computational cost by decreasing the number of stages of the method. 

Depending of the nature of the problem we can try to find the appropriate terms 
in e l r p that must vanish in order to have an optimal performance. For example, 
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if we consider a method such that the coefficients a%,bi satisfy p^o = Px,x = 1 and 
P2,i = P3,i = P4,i = 0, then, 

\U-{A + eB)\ = (D(et 4 + eV), 

but as |e| <C r this method is of effective order 4. In a more general context we 
will have methods such that, 

\H-(A + sB)\ = 0(et Si + eV 2 + eV 3 + ■ • ■ + £ "V m ). (18) 

We remark that si is the error of consistency for the scheme, i.e. is the error 
behaviour in the limit case e — > 0. Nevertheless, in many cases for small step-sizes 
the method can behave as one of order S2- In what follows we will refer to the 
generalised order of a method in terms of the order in powers of e. Hence, we will 
say that a method has order (si, S2) if \H — (A + eB)\ = 0(er Sl + e 2 t S2 ). In terms 
of the local error, we have |5(r) - exp[r(^ + eB)} \ = C(er Sl + 1 + e 2 r S2 + 1 ). 



5.1 ABA schemes of generalised order (2n,2) 



McLachlan (1995) noted that as |e| -C r we can have high-order methods by only 
killing the terms in er k . Independently 



Chambers and Murison 



Robutel (2001 1 dealt with this problem following similar ideas, 



(2000); 



Laskar and 



Laskar and Robu- 



tel 2001 1 providing an explicit computation of the coefficients of the remainder 
for all order k. One of the main advantages of only killing the terms in er k is that 
we are sure that all the coefficients <ij,6j will be positive. As a consequence the 
coefficients aj,6j will be small and the numerical scheme will be stable. 

In Table [2] we summarise the coefficients for the different ABA (2n, 2) schemes 
for n = 1, . . . , 4. For further details on how to find the aj, 6j coefficients and the 
coefficients for n > 4 see ( |McLachlan 1995 Laskar and Robutel[ 2001). Since all 
the methods we consider are symmetric, we only collect the necessary coefficients 
of each scheme. Thus, ABA(8, 2) corresponds to the composition 



gQi-rA gbierS ^etB ^tA ^etB ^itA ^>\etB ^xtA 



We will follow this convention throughout the text. 

In Figure[l]we compare the performance of the ABA(2n, 2) for n = 1, 2, 3, 4 for 
the 4 inner planets (left) and the 4 outer planets (right). The a;-axis corresponds 
to the cost of the scheme (r/s) and the t/-axis corresponds to the maximum energy 
variation for one integration at constant step-size r. iLaskar and Robutel (20011 



already saw that the optimal schemes for this problem where those of orders (6, 2) 
and (8,2) (i.e. SABA3 and SABAi following their notation). 

The error on the Hamiltonian approximation of these schemes is 0(eT 2n +e 2 r 2 ). 
In Figure [l] we can see how the error in energy decreases in r with slope In for 
large steps-sizes and slope 2 for smaller steps-sizes. We also see how for small step- 
sizes there is no difference between the cost of the ABA62 and .4,6.482 schemes. In 
order to improve their performance we need to kill the term in e 2 r 2 rather than 
those of order er 2k for k > 4, which are the limiting factor of these schemes. 
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Table 2 Coefficients for the ABA(2n,2) methods for n = 1,...,4 (Laskar and Robutcl 
[200l| l. 



id 


order 


stages 


ai, bi 




ABA22 


(2,2) 


1 


ai = 
6j = 


1/2 
1 


ABA42 


(4,2) 


2 


ai = 

«2 = 
bi = 


1/2 - x/3/6 
1/2 


ABA62 


(6,2) 


3 


ai = 
ai = 
bi = 
b 2 = 


1/2 - v^/lO 

VT5/10 

5/18 

4/9 



ai = 1/2 - ^525 + 70^/70 

a 2 = ( ^525 + 70^30 - ^525 - 70^) /70 

ABA82 ( 8 ' 2 ) 4 a 3 = ^525 - 70^/35 

hi = 1/4 - v / 30/72 

fe 2 = 1/4 + v / 30/72 



5.2 .46.4 schemes of order (2n,4) 

In this section we will describe three different procedures to cancel the dominant 
term e 2 r 2 in order to get methods of generalized order (2n, 4), and discuss their 
performance for the different test models described in Section [4j 



5.2.1 The corrector term (SC) 



Since in Jacobi coordinates A is quadratic in p and B depends only of q, then it 
follows that the term [[A, B], B] depends only on q and thus exp(r 3 e 2 [[A, B], B]) 
can be easily computed. Laskar and Robutel (2001) noticed that it is possible to 



incorporate this term into the previous compositions with a conveniently chosen 
constant so as to cancel the term of order e 2 r 2 in the asymptotic expansion Eq.|l5[ 
We note that this corrector scheme is different than the one introduced by Wisdom 
et al ( 1996 1 where the corrector added at the beginning and at the end of each 



step-size is a change of variables. 

Thus, l et S n (t) be one of the symplectic ABA schemes of order (2n, 2) described 
We can get rid of the term in e 2 r 2 by considering 



in Section 



5.1 



5C„(r) =exp(-rV^p,B],B]) 5„(r) exp (-rV^p, B], B]) 



(19) 



with the appropriate choice of the constant c. In Table [3] we show the value for 
the coefficient c for each the four ABA{2n, 2) schemes described before. For further 
details see (Laskar and Robutel 2001 1. Notice that SC n corresponds to integrating 
log(5„(r)) 

we obtain 



T^e^cL sia,b\,b\ using the leapfrog scheme. 
So using Eq. 



19 



with any of the _4£>.4 (2n, 2) scheme in Section 



5.1 



a new integrating scheme of order (2n, 4) with no negative intermediate step. 
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Fig. 1 Comparison between the ABA(2n, 2) methods for n = 1, 2, 3, 4 applied to the 4 inner 
planets (left) and the 4 outer planets (right). The x-axis represents the cost (t/s) and the 
?/-axis is the maximum energy variation over one integration with constant step-size t. Here s 
is the number of stages. 

Table 3 Coefficients c for the corrector term applied to the ABA(2n, 2) schemes in Table [2] 
jLaskar and Robutel)|2001| . 

order c 



1 1/12 

2 (2 - \/3)/24 

3 (54-13^)7648 

4 0.003396775048208601331532157783492144 



5.2.2 The composition scheme (S2 m ) 



Yoshida ( 1990 1 and Suzuki ( 1990 1 independently came up with the same idea to 



find a symmetric scheme of order 2k + 2 from one of order 2k. They both noticed 
that if 5(r) is a scheme of order 2k, then: 



S2k{r) = S(xqt)S{x\t)S(xqt) , 



(20) 



is a scheme of order 2k + 2 for an appropriate choice of the constant coefficients 
xo, x\. One can check that xq, x\ must satisfy 2:eo+:ei = 1 and 2x^ j k+1 +x 2 i k+1 = 0. 
Notice that the second condition is used to cancel all the terms of order 2k, while 
the first one is only for consistency. 



Laskar and Robutel (2001 1 used this idea to turn any of the ABA schemes of 



order (2n, 2) into one of order (2n, 4). If S{t) is a symmetric ABA scheme of order 
(2n,2) then the composition: 



52 m (r) = S m (y T)S(yir)S m (y T). 



(21) 
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is a symmetric method of order (2n, 4) if yo, yi satisfy 2myo + y% = 1 and 2mj/Q + 
yl = so, (y = V(2m - (2m) 1 / 3 ), Vl = -(2m) 1 / 3 / (2m - (2m) 1 / 3 )). 

We have done several test to determine the optimal value of m, and our test 
show that this one is given by m = 2. These results are consistent with those 



of (Suzuki 



1990 



McLachlan 2002 1 who did a similar study in a more general 
framework. The main advantage of this scheme is that we can use it for both 
Heliocentric and Jacobi coordinates. 



5.2.3 McLachlan extra stage scheme (ABA84) 



McLachlan ( 1995 ) studied the possibility of adding an extra stage to the ABA(2n, 2) 
schemes to get rid of the e 2 r 2 term. To add an extra stage derives into having an 
extra pair of coefficients a$ , 6j and an extra algebraic equation to satisfy. All the co- 



efficients will no longer be positive (Suzuki 1991 1. In general, if the coefficients are 



not very large, these methods are stable. The coefficients for the ABA method of 
generalised order (8,4) provided by McLachlan (1995) are summarised in Table [4j 



Table 4 Coefficients for the ABA method of order (8, 4) found by ( McLachlan 1995 i 



id 



order stages 









ai = 





07534696026989288841652780368 








a-2 = 





51791685468825678230077397850 


ABA84 


(8,4) 




13 = 


-0 


09326381495814967071730178218 


5 


bi = 
b 2 = 
b 3 = 





-1 


19022593937367661924523076274 
84652407044352625705508054465 
07350001963440575260062261477 



5.2.4 Results 

In Figure [2] we compare the performance of these three different approaches to 
build methods of generalised order (8,4) against the .46.482 scheme. In the plots 
we show the cost (r/s) vs the maximum energy variation for the three test models: 
the 4 inner planets (left), the 4 outer planets (middle) and the 8 planets in the 
Solar System (right). 

As we can see, the three different schemes improve considerably the perfor- 
mance of the .46.482 (red line). In all cases the corrector scheme SC (blue line) 
and the McLachlan .46.484 (purple line) show a similar quantitative behaviour. 
The difference between them is the cost of the extra stage in .46.484, as we are 
assuming that the corrector is completely free. We note that this is not entirely 
true if the number of bodies is large (n > 4). On the other hand, the composition 
methods, Se m (green line), improves the performance with respect to the .46.482 
(red line) but is much more expensive than the other two schemes. 



5.3 ABA schemes with generalised order (si, S2, . . . ) 



In the previous section we have seen that adding an extra stage to cancel the term 
of order e 2 r 2 gives good results. We can extend this idea and add more stages to 
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Fig. 2 Comparison between the different schemes to kill the the e 2 r 2 terms in the .46.482 
scheme. From left to right: the 4 inner planets, the 4 outer planets and the whole Solar System. 
The x-axis represents the cost (t/s) and the y-axis the maximum energy variation for one 
integration with constant step-size t. 



kill the error terms accounting to the main limiting factor for each problem ( Blanes 



et al 2012). This translates into adding more constraints on the coefficients. As 



long as the increase in the computational cost is less than the gain in accuracy 
these methods will be competitive. In Figure [2] we see that the dominant error 
term for the ABA84 varies between the different test models. Notice that for the 
outer planetary system the scheme behaves as one of order 4, so the dominant term 
is £ 2 r 4 . On the other hand, for the inner planetary system, the method behaves 
as one of order 8, now the dominant term is er 8 . 

Hence, to improve the performance of the McLachlan's ABA84 we need to kill 
different error terms depending on the problem. For the inner planets, a method 
of order (f 0, 4) should perform better than one of order (8, 6). While for the outer 
planets a method of order (8, 6) should give better results than one of order (10, 4). 

(20121 we find details on how to solve the algebraic equations 



In 



Blanes et al 



and find the set of coefficients a,,6i that provide ABA schemes for a given arbi- 
trary order (si,S2, - • •)• We must mention that there is no unique combination of 
coefficients dj, b% for a given order. From all the possible solutions we have selected 
those that give a better approximation and whose coefficients aj,b, are small. In 
Table [5] we summarise the coefficients for three ABA schemes: one of order (10,4); 
one of order (8,6,4) and one of order (10,6,4). 



5.3.1 Results 



In Figure [3] we compare the performance of the three schemes summarised in 
Table [5] against the ABA82 and ABA8A for the three test models. 

In the left-hand side of Figures[3]we have the results for the 4 inner planets. We 
recall that the dominant error term in the ABA84 scheme was er 8 . Hence, a method 
of order 10 in e should perform better than one of order 8 in e. Nevertheless, as we 
can see there is no significant gain in the performance of these schemes with respect 
to ABA8A-. Apparently, for these methods the gain in precision is proportional to 
the computational cost in this range of accuracy. 
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Table 5 Coefficients for A BA symmetric splitting methods of orders (10,4), (8,6,4) and 
(10,6,4) planes et al||2012) . 



id order stages 









ai 


— 





047067100645972506129478876372 










= 





184756935417088106924737619370 















282706005679836205324361656554 


ADA 1 A/1 


/in a \ 
(1U,4) 


i 


a 4 


= 


-0 


014530041742896818378578152296 


bi 


= 





118881917368197019945350395085 








b-2 


= 





241050460551501565744166786590 








b 3 


= 


-0 


273286666705323806054311398166 








bi 


= 





826708577571250440729588432981 








ai 


= 





071133426498223117777938730006 








ai 


= 





241153427956640098736487795326 








0.3 







521411761772814789212136078067 


ABA864 


(8,6,4) 


7 


Q>4 




-0 


333698616227678005726562603400 


bi 







183083687472197221961703757166 








b 2 







310782859898574869507522291054 








b 3 




-0 


026564618511958800697212137916 








bi 







065396142282373418455972179391 








ai 







038094497422412195456975322308 








ai 







145298716116913749294020072660 








0.3 







207627695725541250716205611324 








°-i 







435909703651526159223154862401 


ABA 1064 


(10,6,4) 


8 


05 




-0 


653861225832786709380711737390 








bi 







095858880837075210610771503771 








b 2 







204446153142998780680507783916 








b 3 







217070347978991101714338592430 








bi 




-0 


017375381959065093005617880118 



In the middle of Figure [3] we have the results for the 4 outer planets. We recall 
that here the dominant term in the .46.484 scheme was e 2 r 4 , hence we expect the 
schemes of order 6 in e 2 to be better than the .46.484. As we can see the .46.4864 
and the .46.41064 schemes give better results that the .46.484. In both cases the 
optimal cost is around 1CF 2 vs an optimal cost of around 1CP 3 for the .46.484 
scheme. 

The main difference between the inner planets and the outer planets is the size 
of the perturbation. We recall that in Jacobi coordinates, for the inner planets 



e « 4.54 ■ 10~ b , while for the outer planets e w 2.03 • 10" 4 (see Table [lj. The 
difference of about 2 orders of magnitude should explain the difference in the 
performance of the different schemes, as the relevance of the terms eV in the 
error approximation will vary depending on the size of e. 

Taking this into account, one can be surprised by the performance of these 
schemes when we consider the whole Solar System (Figure [3] right). Here the size 
of the perturbation (e £s 1.96 ■ 10~ 4 ) is of the same order of magnitude as the case 
of the outer planets. But as we can see in Figure [3] the schemes behave in the 
same way as the case of the inner planets, where the terms of order er 8 dominate 
those of order e 2 r 4 . We think this is due to Mercury: its fast orbital period and 
large eccentricity is limiting the performance of the methods. This phenomena was 
already noticed by Wisdom et al (1996) and re-discussed by Viswanath (2002 1. 
This is why Saha and Tremaine ( 1994 1 proposed to use independent time-steps for 
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Fig. 3 Comparison between the ABA splitting schemes of arbitrary order (si,S2,S3) sum- 
marised in Table [5] against the .46.482 and .46.484. From left to right: the 4 inner planets, 
the 4 outer planets and the whole Solar System. The x-axis represents the cost (t/s) and the 
j/-axis the maximum energy variation for one integration with constant step-size r. 

each planet. They used the leapfrog scheme and adapted it to take fractions of the 
given step-size for each planet, depending on their orbital period. It is not trivial to 
extend these ideas using the higher order schemes described in this section, and a 
second order method is not the appropriate option to achieve round-off accuracy. 



6 Splitting Symplectic Integrators for Heliocentric 

We recall that all the tests done in Section [5] have been done using Jacobi coordi- 
nates. All these integrating schemes assume that the two parts of the Hamiltonian 
Ha,Hb are integrable. This is true for Jacobi coordinates where: 

H Jb (q,p)=H K (q,p) + Hi{q). (22) 

where Hx(q,p) is integrable (it is a sum of independent Kepler problems) as well as 
Hj(q) (it only depends on q). However, this is not true for Heliocentric coordinates 
where: 

H He (q,p) = H K {q,p) + Hj(q,p), (23) 

and Hj(q,p) is not integrable, which can be a problem if we want to apply the 
splitting schemes presented in Section [5] 

An option to deal with the non-integrability of Hj(q,p) is to use another nu- 
merical method to integrate Hj{q,p) and compute the exp(&iT-B) up to machine 
accuracy. Unfortunately, this can drastically increase the computational cost of 
the schemes. 

We propose to use the fact that Hj(q,p) = Ti(p) + Ui(q) splits naturally into 
two parts, one depending on positions, the other in velocities. These two parts are 
integrable on its own and small with respect to Hk (l, p)- I n & general framework, 
the Hamiltonian splits as: 



H = H A +e{H B +H c ), 
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where Ha, Hg and He, are integrated when we consider them separately. In the 
same way as in Section [5] we could try to find appropriate coefficients a,-, b^, Cj, 
such that 

s 

S(t) = Y\ exp(airA) exp(ebjrB) exp(ecirC), 

approximates exp(rZ/#). As before, to simplify notation A = {Ha, ■}, B = {Hg, ■} 
and C = {Hq,-}- Then one has to deal with the Lie Algebra generated by A, B 
and C. The number of order conditions as well as the complexity to solve them 
numerically to get the coefficients a,i,bi,Ci grows extraordinarily with the order 



Blanes et al (2012). A simple alternative way to proceed is to use the splitting 



ymplectic schemes in Section [5] 

a 

S(t) = Yl exp^TA) exp(eb l r(B + C)) 
»=1 

and approximate exp(ebir(B + C)) with 
exp(ebiT(B + C)) 



exp(e^rC) exp(efo i r_B) exp(e^rC). 



(24) 



(25) 



Here we take C as the Lie operator associated to T\ (p) due to its lower computa- 
tional cost. 

In general Hg and Hq do not commute {{Hb,Hq} ^ 0), so this approximation 
adds an extra error contribution term, e 3 r 2 , which will be negligible for small e. 
In Figure [2] we see the result of taking the ABA&2, the ABA84 and the S2 m 
splitting schemes using Eq. [25] to deal with Heliocentric coordinates. We can see 
that in general the symplectic schemes have the same behaviour as with Jacobi 
coordinates (Figure [2j. 



Mer - Ven - Ear - Mar (Helio Coord) 



Jup - Sat - Ura - Nop (Helio Coord) 



Mercury to Neptune (Helio Coord) 






Fig. 4 Comparison between the ABA&2, .46.484 and S2 m schemes discussed in Section ! 
applied to Heliocentric coordinates. From left to right: the 4 inner planets, the 4 outer planets 
and the whole Solar System. The x-axis represents the cost (t/s) and the y-axis the maximum 
energy variation for one integration with constant step-size t. 



We do see a difference in the case of the outer planets (Figure [4] middle) . Now 
the ABA84 scheme behaves as one of order 2 for small step-sizes. This is due to 
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the extra error term e 3 r 2 . We recall that the main difference between the inner 
and the outer planets is the size of the perturbation (e) which is smaller in the 
first case. Here the terms of order e 3 t are negligible for the inner planets but not 
for the outer planets. 

Unfortunately, when we consider high-order symplectic schemes like the ones 
presented in Sections |5 .3| these extra error term will become relevant, jeopardising 
the performance of these schemes. 

Chambers] ( |1999| ; Wisdom (2006) proposed to rewrite the Hamiltonian in He- 



liocentric variables so that Hg and He commuted ({Hg,Hc} = 0), and then used 
the fact that exp(e6 i r(B + C)) = enp(ebiTB) enp^ebirC) . For further details see 
Appendix [Cl 



6.1 ABAH (2n,4) specific for Heliocentric coordinates 



We have just seen that for Heliocentric coordinates we can adapt the splitting 
schemes described is Section [5] using Eqs. |24| and |25| But with this an extra term 
e t appears in the error approximation that will limit the performance for high- 
order schemes. One can check that this error terms is associated to the algebraic 
expression b\ + b% + ■ ■ ■ + b^ . We can add an extra stage to the scheme so that it 



also satisfies: 



b? + bl + ■ 



+ b" 



0, 



leading to symplectic schemes with the same generalised order as before for He- 
liocentric coordinates. Table [6] collects the coefficients for the ABAH scheme of 
order (8, 4) for Heliocentric coordinates. This scheme has the same effective order 
as the McLachlan ABA scheme of order (8,4) (Section 5.2.3) but it is specific 
for Heliocentric coordinates. In Figure [5] we compare the performance of this new 
scheme against the _4£>.484 scheme. As we can see, for the outer planets the new 
ABAH844 scheme behaves better that the A6.484 for small step-sizes. 



Mer - Ven - Ear - Mar (Helio Coord) 



Jup - Sat - Ura - Nop (Helio Coord) 



Mercury to Neptune (Helio Coord) 




-3 -2 -1 





Fig. 5 Comparison between ABA84 and ABAH844. From left to right: the 4 inner planets, 
the 4 outer planets and the whole Solar System. The x-axis represents the cost (t/s) and the 
j/-axis the maximum energy variation for one integration with constant step-size r. 
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Table 6 Coefficients for ABAH specific sym metric splitting me thods for Heliocentric coor- 
dinates of orders (8,4), (8, 6, 4) and (10, 6, 4) l |Blanes et aj||2012) . 



id 


order 


stages 
















a i 


= 


0. 


,27414026894340187616405654402 








a -2 


= 


-0 


. 10756843844016423062511052968 








0-3 


= 


-0 


. 04801850259060169269119541721 


ABAH844 


(8,4) 


6 


0,4 


= 





, 76289334417472809430449880574 















, o40oob 1 yblb^ol^ ( 1 / Mzz4yilb49 








b 2 




-0 


, 85857544895678285658812832469 








b :i 







. 71768965379427013885587920820 








a\ 




0. 


, 06810235651658372084723976682 








ai 







, 25113603872210332330728295804 








0-3 




-0 


. 07507264957216562516006821767 








a 4 




-0 


.00954471970174500781148821895 


ABAH864 


(8,6,4) 


8 


0-5 







. 53075794807044717763406742353 








In 







, 16844325936189545343103826977 








b 2 







, 42431771737426772243003516574 








b 3 




-0 


, 58581096946817568123090153554 








b 4 







. 49304999273201250536982810002 



. 04731908697653382270404371796 
. 2651 1052357487851595394800361 
. 00997652288381124084326746816 
. 05992919973494155126395247987 
. 25747611206734045344922822646 
. 1 1968846245853220353128642974 
. 37529558553793742504201285376 
. 46845934183259937836508204098 
. 33513973427558970103930989429 
0.27667111912108009750494572633 



ai 
a 2 
a3 

a 4 

ABAH1064 (10,6,4) 9 ^ 

b 2 
63 
64 
br, 



6.2 ABAH specific methods with arbitrary order (si, s2, ...) 



In the same way we can add the extra constraint b\ + ■ ■ ■ + bn = to the high-order 
schemes in Section [573| to have high order splitting schemes specific for Heliocentric 
coordinates. In Table |6] we show the coefficients of two ABAH schemes of orders 
(8,6,4) and (10,6,4). All these schemes have one more stage than the schemes 
presented in Table [5] 

In Figure [6] we compare the performance of the .4,6.482 with the other three 
schemes in Table |6j Where the behaviour of the schemes depending on its order is 
similar to the one presented in Jacobi coordinates. For the inner planets (Figure]^ 
left) all ABAH schemes present a similar optimal cost. For the outer planets (Fig- 
ure [6] middle) the ABAH schemes of order (8,6,4) and (10,6,4) are much better 
than the other two schemes. We recall that here the size of the perturbation is 
larger and killing the terms of order e 3 r 4 does make a difference. Finally, if we 
consider the 8 planets in the Solar System (Figure|6]right) here the ABAH schemes 
of order (8, 6, 4) and (10, 6, 4) do improve the performance of the schemes of order 
(8,4). We recall that this was not the case in Jacobi coordinates (Figure[3|. 
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Fig. 6 Comparison between ABJKH schemes of order (8,4,4), (8,6,4) and (10,6,4) specific 
for Heliocentric coordinates and the .4B.482. From left to right: the 4 inner planets, the 4 outer 
planets and the whole Solar System. The tr-axis represents the cost (r/n) and the j/-axis the 
maximum energy variation for one integration at constant step-size r. 



7 Jacobi vs Heliocentric coordinates 

In Sections [5] and [6] we have described different splitting schemes for both Jacobi 
and Heliocentric set of coordinates. In the case of Heliocentric coordinates the 
expressions for the Hamiltonian are less cumbersome and easier to handle (see 
Appendix [b]). But the size of the perturbation is larger than in Jacobi coordinates 
and an extra stage to deal with the non-integrability of Hj must be added to have 
high-order schemes. Here we want to compare the performance of the different 
schemes for both set of coordinates. 

We compare the performance of the ABA methods of order (8,2), (8,4) and 
(10, 6, 4) in both set of coordinates, with the three test models used throughout 
this report. We recall that the (8,4) and (10,6,4) schemes have an extra stage 
in Heliocentric coordinates. In Figure [7] we summarise the performance of these 
schemes. From left to right we have the results for the inner planets, the outer 
planets and the whole Solar System. We distinguish the order of the schemes 
by the colour. Where the lines in red represent the schemes of order (8,2), the 
blue lines those of order (8,4) and the purple lines those of order (10,6,4). We 
use continuous lines to refer to Jacobi coordinates and discontinuous lines for 
Heliocentric coordinates. 

If we look at the results for the inner planets (Figure[7|left), we can see there is 
not much difference between taking Jacobi or Heliocentric coordinates (continuous 
vs discontinuous lines). In both cases the size of the perturbation is small (Table[lJ 
and there is not much difference between taking a splitting scheme of order (8, 4) 
or (10,6,4). In both cases the terms in e in the error expansion are the ones that 
matter, but there is not much difference between taking order 8 or 10 in er k . 
We should have to use arithmetics with higher precision to see the difference (see 
Appendix [d]) . Hence the ABA84: is the best choice for this case. 

If we look at the results for the outer planets (Figure [7] middle) , again there 
is no significant difference between Jacobi and Heliocentric coordinates. But here 
the ABA schemes of order (10, 6, 4) performs much better that the other schemes, 
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Mercury - Verus - Earth - Mars Jupiter - Saturn - Uranus - Neptune Mercury to Neptune 




-5 -4 -3 -2 -1 -5 -4 -3 -2 -1 -5 -4 -3 -2 -1 



Fig. 7 Comparison between Jacobi (continuous lines) and Heliocentric (discontinous lines) 
coordinates using the schemes ABA82 (red), .46.484 (blue) and .40.41064 (purple). From left 
to right: the 4 inner planets, the 4 outer planets and the whole Solar System. The x-axis 
represents the cost (r/s) of the method and the y-axis the maximum energy variation for one 
integration with constant step-size r. 

having an optimal step-sizes one order of magnitude larger than the ones for the 
schemes of order (8,4). 

If we look at the whole Solar System (Figure [7] right), we see that there is a 
big difference between taking Jacobi or Heliocentric coordinates. Looking at the 
ABA82 scheme (red lines) we see that the slopes are the same but that there is a 
difference of about one order of magnitude in accuracy for a given step-sizes. If we 
look at the scheme of order (8,4) (blue lines), we see that in Jacobi coordinates 
the methods behaves as one of order 8, while in Heliocentric coordinates this one 
behaves as one of order 4. This can be explained by the difference in the size of 
the perturbation (see Table [T]) in both set of coordinates. We also see that there is 
a big difference between the optimal step-size for both set of coordinates, making 
Jacobi coordinates by far the best choice. Finally, if we compare the .46.4 schemes 
of order (10,6,4) (purple lines) the difference between the two set of coordinates 
is drastically reduced, although Jacobi coordinates still perform slightly better. 
While the extra stages to go from an (8,4) scheme to a (10,6,4) one do not 
improve in Jacobi coordinates. This is not the case of Heliocentric coordinates, 
where the (10,6,4) gives the best results and the difference between the two set 
of coordinates is not as relevant as before. 

Although Jacobi coordinates presents better results for most of the test models, 
we believe that using Jacobi or Heliocentric coordinates is a matter of choice. 



8 Conclusions 

In this article we have reviewed different symplectic splitting schemes and tested 
their performance for the case of the planetary motion. We recall that in the case of 
the planetary motion, using an appropriate change of variables, the Hamiltonian 
of the N - body problem can be rewritten as Hr + Hj. A sum of independent 
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Keplerian motions for each planet (H^) and a small perturbation term given by 
the interaction between the planets (Hi). 

There are two set of canonical coordinates that allow us to write the Hamil- 
tonian in this way: Jacobi and Heliocentric coordinates (Section [3j, Although in 
Jacobi coordinates the size of the perturbation is smaller and Hj is integrable, 
Heliocentric coordinates seem more natural and the expressions are easier to han- 
dle (Appendix [B]) . In this article we have compared the performance of different 
symplectic splitting schemes in both set of coordinates. 

In Section [5] we described different splitting symplectic schemes for Jacobi 
coordinates. In Section [6] we saw how to extend these schemes to use Heliocentric 
coordinates. We note that all the splitting schemes for Jacobi coordinates can also 
be used in Heliocentric coordinates, but in order to have a comparableperformance 
an extra stage to kill the terms in e 3 t 2 must be added (see Section ml). 

We have seen that in Jacobi coordinates, the ABA84. scheme introduced by 



McLachlan (19951 and the ABA1064 scheme Blanes et al (2012) give the best 



results when we look at the motion of the whole Solar System. The high eccentricity 
of Mercury and its fast orbital period are the main limiting factors and taking 
higher order splitting schemes do not always provide significant improvments. But 
for different planetary configurations, as the 4 outer planets, the ABAW6A has a 
better performance than the ABA8A. 

When we consider Heliocentric coordinates, the ABAHW64: (Blanes et al 2012 1 
gives the best results when we consider the whole Solar System. In this case, 
probably because the size of the perturbation is larger, adding extra stages to 
have higher order schemes does improve the results. 

Moreover, the performances of the schemes in both set of coordinates, Jacobi 
or Heliocentric are very similar for the scheme of order (10,6,4), with a slight 
advantage for the Jacobi coordinates. Depending on the problem, one can thus 
use either system of coordinates, but it is clear that using high order schemes as 
the ABA(n)864: and A8.4('H)1064 |Blanes et al| ( |2012[ ) can drastically improve the 
results. This should be even more the case for highly perturbed systems as some 
extra solar planetary systems with close planets of large masses. 
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A On the Compensated Summation 

Using any of the symplectic integrating schemes described in this article, we require successive 
evaluations of exp(a^rA) and cxp(birB). Each of these evaluations slightly modifies the posi- 
tion and velocity of each planet. For r small, we will have a loss in accuracy due to round-off 
errors. The Compensated, Summation is a simple trick that is commonly used to reduce the 
round-off error. In a general framework, when we consider a numerical method for solving an 
ODE, we require a recursive evaluation of the form: 

y n+1 = y n +8 n , (26) 

where y n is the approximated solution and 8„ is the increment to be done. Usually <5 n will be 
smaller in magnitude than y n . In this situation, the rounding errors caused by the computation 
of 8 n are in general smaller that those to evaluate Eq. |26| The algorithm tha t can be used in 
order to reduce this round-off error is called the "Compensated Summation" | |Kahan[ |1965} . 

Compensated Summation Algorithm: Let yo and {<5 n } n >o be given and put e = 0. 
Compute 2/1,2/2, •■■ from Eq. \ 26\ as follows: 

for n = 0, 1,2, . .. do 

a = y n 

e = e + <5 n 

y-a+l = a + e 

e = e + (a - y n +i) 
enddo 



This algorithm accumulates the rounding errors in e and feeds them back into the sum- 
mation when possible. At each time-step of the integration, when we evaluate exp(airA) or 
exp(feiT_B), the increment in position and velocity is done using the compensated summation. 
In Figure [s] we show the results for the ABA (2n,2) schemes for n = 1,2,3,4 using double 
(left) and extended precision (right). In both cases we gain almost one order of magnitude in 
precision when we take into account the compensated summation. 



Jup-Sat [Double Precision CS Vs NO CS] Jup-Sat [Extended Precision CS vs NO CS] 




-6 -5 -4 -3 -2 -6 -5 -4 -3 -2 



Fig. 8 Comparison between the ABA schemes of order (2n, 2) for n = 1,2,3,4 applied to 
the Sun-Jupiter-Saturn three body problem. With (CS) and without (noCS) the compensated 
summation. The x-axis represent the cost (r/n) and the j/-axis the maximum energy variation 
for one integration with constant step-size r. Left: using a double precision arithmetics; Right: 
using an extended precision arithmetics. 
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B Integration Schemes (some help on the practical coding) 

In this paper we have reviewed many splitting symplectic integrating schemes, all of them of 
the form: 

5(t) = exp(airj4) exp(bir_B) ... exp(bi tB) exp(airA), (27) 

where exp(rA) and exp(rB) can be computed explicitly. They correspond to the integrals of 
the two different parts of the original Hamiltonian. In this section we show how to compute 
explicitly exp(Tj4) and exp(rS) for the particular case of the N-body problem in Jacobi and 
Heliocentric coordinates. We note that from now on: u stands for the momenta associated to 
u and u' stands for du/dt. 



B.l Keplerian Motion (H K ) 

We recall that in Jacobi coordinates, 



v 2 m- i 

whereas in Heliocentric coordinates, 



»K = J2(l — ^ L -G m ^), (28) 

i=l 





" ' m + mi 


(>[[' 






_ m mi _ 



G Tr-rr • ( 29 ) 



i=l 

In both cases Hk is a sum of independent Keplerian motions. In Jacobi coordinates each 
planet follows an elliptical orbit around the centre on mass of the Sun and the planets that 
are closer to the Sun, the mass parameter of the system is fij = Gm. While in Heliocentric 
coordinates each planet follows an elliptical orbit around the planet-Sun centre of mass, and 
the mass parameter of the system is Hh = G(mo + mi). 

It is well known that Kepler problem is integrable, but the solution from time t = to to 
t = to + t is expressed in a simple form if we consider action-angle variables. To compute 
exp(rL[j K ) we need to be able to compute (r(trj + t), v(trj + r)) from (r(to), v(to)). 

An option is to change to elliptical coordinates, advance the mean anomaly and then 
return to cartesian coordinates. But this can accumulate a lot of numerical errors as well as it 
is very expensive in terms of co mputational cost. Instead we use a similar idea as the Gauss / 
and g functions (Danby 1992), where we use an expression for the increment in position and 
velocities for a given step-size t, without having to perform any change of coordinates. Let us 
give some details on how to derive these expressions. 

In elliptical coordinates the motion of the two body problem is given by (a, e, i, Q, u, E), 
where all of the elements remain fixed except for E that varies following Kepler equation 
(n(t — t p ) = M = E — es'mE). Using a reference frame where the orbital plane is given by 
Z = 0, the Jf -axis is the direction of the perihelion and the y-axis completes an orthogonal 
reference system on the orbital plane, the position (X, Y, 0) and velocity (X' ,Y' , 0) are given 
by: 

X = a(cos E — e), Y = a\/l — e 2 sin E, 

2 2 ( 3 °) 

X> = sin£, Y' = — A /l^ e 2cosE, 

r r 

where r = a(l — ecosi?) and n = /x 1 / 2 a — 3 / 2 . The position and velocities on a fixed reference 
frame are given by: 

«'\ I XX' 

yy'\= n 3 (n) x ni{i) x rc 3 M x [y y> I , (31) 



where 







10 \ / cos sin 

7^(0)= | cos6» sin6» , and 7e 3 (6») = -sin6» cos6» 
- sin 6 cos 6 / \ 1 



High precision Symplectic Integrators for the Solar System 



25 



we have that: 



Notice that 

K 3 (f2) x Ki(i) x K 3 (u) = Kx K 3 (m), 
where m = Q + u> and Tl = K^tt) x Ki(i) x Tl 3 {-Q). Given that = 7li(i/2)7l 1 (i/2) 

1 - 2p 2 2pq 2 PX 
2pq 1 - 2q 2 -2q X 
-2p X 2q X 1 - 2p 2 - 2q 2 ; 

where p = sin i/2 sin Q, q = sini/2sini?, and \ = ~J\ — p 2 — q 2 = cosi/2. From Eq. 1 3 1 1 we 
have, 



7v : 



(32) 



Hence, 



[ r(t ), v(t ) ] = H x 7e 3 (ro) x 
r(t + 5t), v(t +(5t) ] = Tl X 7e 3 (ro) x 



r(t + <5t), v(t + <5t) ] = [ r(t ),v(t ) ] 
= [r(t ),v(t ) ] 



X X' 


Yl Y[ 




Xq Xq 


-l 


~Xi X[' 


Yo Y{_ 




.Yl Y{ 



an a 12 

0,21 d'22 



(33) 
(34) 

(35) 
(36) 



One can check that, 



an = 1 + (cos(B a - E a ) - 1)- 



ro 



a 3 / 2 

a.21 = — rrr sin(Bi — _Eo) — esini?i + esini?o, 



(37) 



a.12 



r ri 



sin(Bi -Eq), 



0,22 = 1 + (cos(Ei - £ ) - 1) — , 

ri 

where = a(l — ecos Ei) for i = 0, 1. We use Kepler's equation to compute SE = E\ — Eq 
from St = ii — <o- Taking = n(t, — t p ) for i = 0,1, we have that 8E is the solution of 



x — e cos E sin a; — e sin E cos x + e sin E — nSt = 0. 



(38) 



Calling C = cos SE, S = sin SE and c e = ecosEo, se = sinBo we have that ri 
a(l — ce ■ C + se ■ S). Now we can rewrite Eq. |37| as: 

an = 1 + (C-1) — , 
ro 



3/2 

a 21 = 5t + (S~SE)^, 



S 



(39) 



ai2 



a22 = 1 + 



r 0\/a(l — ce ■ C + se ■ S) 
C- 1 



1 - ce ■ C + se • S 

To summarise, given r = r(trj), v = v(to) and defining ro = ||r|| and vq = ||v||. We find 

a = ro/(2 — roii 2 ), ce = e cos Eq = rovg — 1, se = e sin Eq = (r, v)/^/Jia. 

Then we take St and we use Eq. |38| to find SE. Finally we use Eqs. |36| and |39| to find 
r(t +<St),v(to +St). 



26 



Ariadna Farres et al. 



B.2 Jacobi Coordinates 



We recall that in this set of coordinates the perturbation part is given by: 



Hi = U\ = G 



^^viiviH Hnii; o^^lln-rjli; 



(40) 



B.2.1 Computing exp(L^f 7 ): 

U\ depends only on the position, hence the equations of motion are given by, 

dUi 



d dUi 

Vk : 



dt 



9 v k 



dt a v k 



t t . - rn-imi . 

Using v. = Vi we have 



VkM = V k (T ), 



VkW = Vk(ro) - t 



dUi 



As the expressions for dUi/dv k can be a little cumbersome, we compute them separately. 
When wc derive Hj with respect to Vk we must derive 3 main expressions: l/||vi||, l/||r;|| 
and l/||r; — rj|| for i < j. Wc first give the derivatives of these factors with respect to Vk and 
then we will deduce dHi /c9vk for k = 1, . . . , n. 



d 
dv k 

d 
dv k 



(llvill) 

(m) 



d 

<9v k \ i 



^..||3 



r.l|3 



rjll 3 




■tpi,j,k, where Vij.fc 



— 1 if i < j = k, 
else (k < i < j, i < j < k). 



To compute dU\/dv k for k = 1, . . . , n, we consider separately the cases k = 1 and k > 1: 



<9[/i _ ^momi 
d vi 771 



E m »TTTTT3 + E m M 



dv k 



Gm k 



Vk-l 



-Vk-i-. 



Vk 

|v k || 



- mo; 



rk 



|r k | 



mo 



E ■ 

i=fe+i 



r k-rj 

aj ||rk-rj||3 



fc-i 



E • 

i=k+l 
i - rk 



E E • 

i=l j=fc+l 



71 T7? 
r: — r: 3 
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B.2.2 Computing the Corrector: exp(L^ j 4 B j B y) 

In Section |5.2,l| we described a splitting symplectic schemed where a corrector term was added 
at the beginning and at the end of each step-size. The corrector term is given by, 

cxp(-rV^L c ), 

with Lc = ^{{A,B},s} an d c a constant coefficient that depends on the order of the ABA 
scheme. 

In Jacobi coordinates A is quadratic in p and B only depends on q so {{A, B},B} only 
depends on q and {{A, B},B} is integrable. We recall that A = Hjiep = To + t/o and B = 
H pert = U\. Hence, 

{{To + C/o, Ui}, Ui} = {{To, Ui}, Ui}. 

Given that T = V — ilYilL ) wc have, 

f-J r)i-i m i 2 

{2bf£ r l} = f : _2i_^ 1 



{{T ,tfi},tfi} = J])- 



Then the equations of motion for Lc are given by: 
VkM = Vk(TTj), 

v k (r) = Vk(to) + r > 27;-— — — — , 
ovi ovjOVk 

where 71. = — 2* . As before, using v; = — — - — -vi we have 

v k (r) = v k ( t0 ) + r 7fe g 2 — j 

Again the expression for — are a little cumbersome and we first show how to derive the 

<9v ; <9v k 



different parts in — — : v;/| |v;| | J , rj/[[rj[[ 3 and r j — r j / 1 1 rj — r j , 



dUl : vs v. ' '•• :i ..».-) ... - .■• ... - ,-. 
<9v k 

(h,k) (n,h)(n,k) N , 

II 113 ll_ I IK J Si,fc) 



<h,k) (r,-rj,h)(ri-rj,k) ' 
ki-rjIP || ri - rj ||5 



From now on we call Acc(i) = -ji , and 

A -Acc(c)( {h > k) 3 < v " h >< v " k > 



<h,k) Q (ri,h)(ri,k) 



e», s = Acc(s) ( ^-i - 3 



= Acc( S ) ( - 3 <r '-^ h)<ri -^' k) 
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Wc can now give the expressions for Vj, k 



d 2 U± „"io m i 



d 2 U! 

9vi<9v k 

d 2 U! 

<9v k dv k 



G- 



G 



3=2 
mom\m k 



1 " 



Gm fe 



-Vk-iA-k +moQk,k H — 2^ m » e »,fc H 2~~ 2^ m i W k,i,k 

Ik i=k+l r> k i=k+l 



k — 1 k — 1 n 

- ^ m^ i>k>k + —j- ^ ^ m.imj^ iijtk 

i=l ^k j=l j=i + l 



G 



m k mi 



Vk 
fc-i 



1 n 

mo€>l,s — Vk-l&k,l,s H ^ m i( m 0@i,s ~ Vk-A,i,s) 



X^ m *^M,sH m i m j^i 

i=l ^ i=l 7=!+l 



B.3 Heliocentric Coordinates 



We recall that in this set of coordinates the perturbation part is given by: 

f i • f ; _ ^-^ miirii 



H I =T 1+ U 1 = £ — — — - G ^ 



0<i<j<n 



(41) 



_B.3._Z Computing exp(rLr 1 ).- 



Notice that Ti depends only on the momenta (f). Hence, the equations of motion are given 

by, 



d 9Ti 

— r k = 

ctt 9 r k 



™„ 



d 9Ti 

— f k = = 0. 

df. 9 r k 



Finally, 



rk(r) = r k (ro) + r 



mo 



mo 



rk(T) = fk(-ro). 



Computing exp(rL Upert ): 

Notice that U\ depends only on the positions (r). Hence, the equations of motion are given by, 

d. dUi n 
dt d r k 



— i"k 

dt 



dU! 
d r k 



fe-i 



Q E ^ (rk _ rj) _ E ^ (rj _ rk) 

> --1 ^kj j=fc+i ifc , 
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Given that r k = mfcr k , we have: 



rk(r) = r k (r ), f k (r) = r k (T ) - r G 



^ A?. 
vj=i fc J 



(r k 



(r k 



C Heliocentric Coordinates (Alternatives for the set of equations) 

The canonical Heliocentric coordinates used in Section [6] are canonical and the position of 
each body is taken with respect to the position of the Sun. The position and their associated 
momenta are given by: 



uo 



+ u n 



The main difference between Jacobi and Heliocentric coordinates is that in the second set of 
coordinates the kinetic energy is not diagonal in the momenta. Instead we have: 



± n u iii 

2 h m ° 



y- ll r ll 



i ll £2=1 fill" 

2 mo 



(42) 



which can be rewritten as 
1 



E ll»i 



i=0 



1 1 

m rrii 



E - 



0<i<j 



"It) 



(43) 



The extra term due to the momenta of the Sun is added to the perturbation part and 
makes it depend on bot h position and veloc i ties. In Section [3] we used Eq|43| to derive the 
Hamiltonian expression. |Duncan et~al| ( |1998^ |Chambers| ( |1099[ ) ; |Wisdom| \2006\ used Eq. [42] 
instead. Here we discuss the main differences between the two sets of equations and compare 
the performance of the integrators presented in this paper for both expressions. 



C.l Two different expressions for Heliocentric coordinates 

As we know in Heliocentric coordinates the Hamiltonian for an n-planetary system takes the 
form: 

H = H K + T 1 + Ui, 

a sum of Keplerian parts, a quadratic term in the momenta and the gravitational interaction 
between the other planets. Using Eq. |43|we have, 



Ti 
Ui 



m + mi 
m mi 



G 



m m,i 



"'() 



E - 

0<i<j<n 

-a E 



An 

0<i<j<n lJ 



(44) 
(45) 
(46) 



The main advantage of this way to split the equations is that Kepler's third law is satisfied 
for the individual planets: n 2 a 3 = G(mo + rrii). But Hj = Ti + Ui is not integrable and 
{T u Ui}^0. 
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Using Eq. |42| we have the splitting introduced by Chambers (1999), 

T , = i IIEIUfilP^ (48) 

2 mo 

x — . mini; 

u; = -g E -^r- ( 49 ) 

0<i<j'<n lJ 

With this way to split the equations the mass paramete r for th e Keplerian orbits is ii = Gmo for 
all of the planets. On the other hand, T* and Iff (Eqs, |48]49[ ) commute, i.e. {T*, U*} = and 
this is a advantage when we build high order splitting schemes. For simplicity let us consider 
Hb = T\ and Hq = U\ for both expressions. We recall that in Helioce ntric coordinates we 
need to integrate exp(r(i? + C)). Using Chambers' splitting (Eqs. [48)49[ l we have: 

exp(r(B + C)) = exp(rB) exp(rC), (50) 

which can be computed exactly and does not introduce any extra error terms to th e split ting 
schemes discussed in Section [5] Instead using the first splitting expressions (Eqs. |45|46| wv 
used 

exp(r(B + O) « exp( |c) exp(ri3) exp(^C), (51) 

and introduced error terms of order e 3 r 2 . To deal with this in Section [fi] derived splitting 
schemes where an extra stage was added to get rid of these extra error terms. 



C.2 Comparisons between the expressions 

We recall that when we use Chambers expression (Eqs. |47|49| l we use the splitting schemes 
discussed in Section \E\ with 

71 

S(t) = [I ex P( a i TL H K » ) exp(biTL T » ) exp(fe 4 rL [7 » ). (52) 
i=l 

While when we use the classical expression (Eqs. |44|46] l we use the splitting schemes discussed 
in Section l6l with 

n 

■S(t) = Y\exp(a i TL HKiip )exp(b i -L Tl )cxp(b i TL Ul )e^p(bi-L Tl ). (53) 



We compare the ABA schemes of orders (8,2), (8,4) and (10,6,4) for both splitting ex- 
press ions. W e recall that the schemes of order (8, 4) and (10, 6, 4) that use the classical s plittin g 
(Eqs. [44)46 1 have one more stage than the schemes used with Chambers splitting (Eqs. [47)49[ l. 



Figure 9] summarise the performance of the different integrating schemes presented ir 
Sections [5] and [6] From left to right we have the results for the inner planets, the outer planets 
and the whole Solar System. The red lines show the performance of the ABAS2 scheme, the 
green lines are for the .4B.484 schemes and the blue lines are for the ABA1064 schemes. We use 
continuous lines when we consider the classical splitting and discontinuous lines for Chambers 
splitting. 

As we can see, there is no significant difference between using one splitting or the other. 
In some cases one is better that the other. The main advantage that the splitting introduced 
by Chambers is that we do not require and extra stage for a high-order scheme. 
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Fig. 9 Com parison between the two expressions for Heliocentric coordi nates: t he classical ex- 
pression (Eqs. [44]46[ l continuous lines and the Chambers expression (Eq. |47|49| l discontinuous 
lines. For the schemes ABA82 (red), ABA8A (green) and .46.41064 (blue) . From left to right: 
the 4 inner planets, the 4 outer planets and the whole Solar System. The x-axis represents 
the cost (t/s) of the method and the y-axis the maximum energy variation for one integration 
with constant step-size r. 



D Comparison in Quadruple Precision 

As we have discussed throughout the article, in many cases we have seen that despite taking 
higher order methods no significant improvement on the performance of the schemes was 
observed. This is the case of the 4 inner planets in the Solar System, where the size of the 
perturbation is so small that the extra stages to increase the order of the schemes are useless. 
Here the round-off error dominates the terms in e 2 and e 4 . Similar results are also observed 
when we consider the whole Solar System. In order to see an improvement we need to use 
higher precision arithmetics. Here we have repeated the test from Sections [5] and [6] for the 
different integrating schemes using quadruple precision arithmetics. We want to illustrate that 
the different schemes of orders (8,6,4) and (10,6,4) perform better that those of order (8,4). 

In Figures [To] and [Tl] we show the results for the same test models used throughout the 
article for Jacobi and Heliocentric c oord inates respectively using quadruple precision arith- 
metics. For Jacobi coordinates (Figure[lO]l we compare the .46.482, .46.484, .46.4104, ABA864 
and .46.41064 schemes. For Heliocentric coordinates (Figure [TTJ we compare the ABAH82, 
ABAH84, .46.4^864 and A BAH1064. As we c an see in Figure [To] for Jacobi coordinates, 
the .40.4864 and A8.41064 l |Blanes et af| |2012[ l do improve the performance of the .46.484 
(McLachlan 19951. Notice also that for the 4 inner planets (Figurc [To] lcft) and the whole Solar 
System (Figure |10| right) the improvement is achieved for small step-sizes, where the energy 
variation is bellow the machines epsilon for extended arithmetics precision. In Figurc [Tl"1 similar 
results are observed for Heliocentric coordinates. 

From these experiments we see how the .46.4 splitting meth ods of orders (8, 6 , 4) and 
(10, 6, 4) for both set of coordinates improve the performance of the^McLachlan ( 1995 ) .4B.484 
and the |Laskar and Robutei] | |2001^ A8.482. 
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Fig. 10 Comparison using Jacobi coordinates between the .4B.482, AB.484, .4,6.4104, 
.46.4864 and .4B.41064 schemes using quadruple precision arithmetics. From left to right: 
the 4 inner planets, the 4 outer planets and the whole Solar System. The x-axis represents 
the cost (t/s) and the y-axis the maximum energy variation for one integration with constant 
step-size t. 



Mer - Ven - Ear 


- Mar (Hello Coord) 


ABAH82 




ABAH (8,4,4) 




ABAH (8,6,4) 




ABAH (10,6,4) 









Jup - Sat - Ura - Nep (Helio Coord) 

' ABAH82 -^i — ' 
ABAH (8,4,4} — *— 
ABAH (8,6,4) — ■— 
ABAH (10,6,4) — H — 



Mercury to Neptune (Helio Coord) 





-6 -5 -4 -3 -2 -1 



-5 -4 -3 -2 -1 -6 -5 -4 -3 -2 -1 



Fig. 11 Comparison using Heliocentric coordinates between the ABA'HS2, ABA1-i84:, 
ABA1-i864 and ABAH1064 schemes using quadruple precision arithmetics. From left to right: 
the 4 inner planets, the 4 outer planets and the whole Solar System. The z-axis represents 
the cost (t/s) and the y-axis the maximum energy variation for one integration with constant 
step-size t. 
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